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SUMMARY 

This paper describes an accurate economical 
method for generating approximations to the kernel 
of the integral equation relating unsteady 
pressure to normal wash in nonplanar flow. The 
method is capable of generating approximations of 
arbitrary accuracy. It is based on approximating 
the algebraic part of the non elementary integrals 
in the kernel by exponential functions and 
then integrating termwise. The exponent spacing 
in the approximation is a geometric sequence. The 
coefficients and exponent multiplier of the expo- 
nential approximation are computed by least 
squares so the method is completely automated. 
Exponential approximates generated in this manner 
are two orders of magnitude more accurate than the 
exponential approximation that is currently most 
often used for this purpose. Coefficients for 8, 
12, 24, and 72 term approximations are tabulated 
in the report. Also, since the method is 
automated, it can be used to generate approxi- 
mations to attain any desired trade-off between 
accuracy and computing cost. 

NOMENCLATURE 

a|< coefficient in exponential 

approximation 

b exponent multiplier in exponential 

approximation 

bk exponent in exponential 

approximation 

C]k element of matrix of least squares 

normal equations 

d-) elements of right-hand-side vector 

of normal equations 

E(b) least squares penalty function 

E (t ) error in f(t), g(t) - f(t) 

E(s,r) error in F(s,r), 

«° -irt 

/ e [g(t) - f(t)] dt 
s 

f(t) function to be approximated, 

equation (12) 

F(s,r) Computational form of incomplete 

modified Struve function occuring in 
planar part of kernel 

F(s,r) Incomplete modified Struve function 

occuring in planar part of kernel 


g(t) an approximation to f(t) 

G(s,r) Computational form of incomplete 

modified Struve function occuring in 
nonplanar part of kernel 

^(s.r) Incomplete modified Struve function 

occuring in nonplanar part of kernel 

Pk exponent coefficient in equation 

(A2). In this paper either Pk = k 
or pk = 2*/ m . 

r Frequency term occuring in argument 

of exponential factor (r = (u>/U) 

• v'Cyo 2 + z o 2 ) times the 
variable of integration of equations 
(6) and (7)) 

s Lower limit of integration 

(s = si//(y 0 2 + z 0 2)) 

t variable of integration in 

definition of F(s,r) and G(s,r) 


The symbols occuring in equations (1) thru (7) not 
defined above are defined in the references. They 
are not used elsewhere in the paper. Equations 
(1) thru (7) merely illustrate a use for F(s,r) 
and G(s,r). 

INTRODUCTION 

This paper describes a method for computing two 
incomplete modified Struve functions that occur in 
unsteady aerodynamics. These two functions are 
the nonelementary part of the unsteady kernel of 
the integral equation relating lift to downwash in 
three dimensional potential flow. 

The incomplete Struve functions approximated in 
this report occur whenever one attempts to compute 
the three dimensional velocity field induced by 
unsteady doublet distributions. The most common 
occurance is when evaluating the acceleration 
potential kernel of the downwash integral equation 
as in reference 1. However, the same Struve 
functions occur even if the velocity potential is 
used, although the kernel is very different. For 
example, a 24 term approximation developed herein 
is used to compute the wake contribution in a 
velocity potential program called SOUSSA 
(reference 2). 

The method developed in this paper is the 
familiar technique of replacing the algebraic 
factor in the defining integral by an exponential 
approximation and then evaluating the integral in 
closed form. The new exponential approximations 
described in this paper have coefficients computed 
by an easily automated least squares process and 
hence can be used to approximate the incomplete 
Struve functions to any desired accuracy. 
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The paper explains the reason for using 
exponential approximations instead of series 
expansions or direct numerical Integration. It 
also surveys several previously published 
approximations and assesses their accuracy. 

Coefficients for four new approximations are 
tabulated in the paper. The least squares process 
is explained in sufficient detail to enable the 
reader to generate his own approximations to fit 
the accuracy or computing time constraints of his 
problem. 


STATEMENT OF PROBLEM 


Consider the integral equation relating 
unsteady pressure to normal wash in three 
dimensional potential flow 


w(x,s) = 

N 

1/(8tt) ( l //- K ( X , £ ; S , a; w/U,M)p(£,a)d£da) (1) 
n = l J n 


A convenient representation of the kernel, K, 
above is the one given by Harder and Rodden 
(reference 3) 


K = expt-iwXg/UHl^Tj/r 2 + K^/r 4 ) (2) 


Mr 4 exp(-1a)S»/U) 
‘ 2 ’ i - f 2 . ^2,3/2 


B 2 (r 2 + s 2 ) Ms 

■i ; 


R 


TT 


j2_3 


-i 


M r (wr/U) exp ( i ws 2 / U ) 


~RV ♦ s 2 W 


Mr 4 exp(-ia>s 1 /U) 
* R(r Z . ,() 


B 2 (r 2 + s 2 ) Ms. 

+ "R“ 


[2 ? 

R 2 


M 2 r 3 (ur/U) exp (-ius^U) 


+i 


R 2 (r 2 + s 2 ) 1/2 


+ 3 A C 4 exPjl - £S/U ) ds 

Sj (r 2 + s 2 ) 5/2 


( 7 ) 


# 


where x 0 is the distance in the free stream 
direction between the sending and receiving points 
and 


T. = cos (y r - Y$ ) (3) 


T 2 = (z 0 cosy r - y 0 sinY r )(z Q cosy s - y Q sinY s ) (4) 
r = (y 0 2 + z 0 2)1/2 (5) 


where yo and zg are the distance from sending 
to receiving points measured normal to the 
freestream and Yr> Ys are dihedral angles. 


The unit function, u(x), in equation (6) is one 
if x > 0 and is zero is x < 0. The only non 
elementary terms in the above expressions for 
Ki and K 2 are the integrals. They can be 
rescaled to 


F(s,r) = / e _irt (l + t 2 )' 3/2 dt (8) 

s 

and 


G(s,r) = / e' irt (l + t 2 )“ 5/2 dt (9) 

s 


Mr 2 exp(-iois,/U) 

V = - u <V Br) Hr f 7 J ' c 2,1/2 


(r‘ + 


exp{ -iu>s, /U) s„ 2 r/lA 

+ ~2 TUri + / - ■ 2 - ( ~ - 2 -' 3/l ds J W 

(r 2 + s. 2 ) i/2 s. (r 2 + s 2 ) 3/2 


These integrals are called incomplete modified 
Struve functions. In subsonic flow the upper 
limit of integration is infinity as in equations 
(8) and (9) above. In supersonic flow the upper 
limit is variable so the integrals in equations 
(6) and (7) are expressed as differences of 
incomplete Struve functions. The algebraic part 
of equations (8) or (9) could be approximated by 
an exponential function. However, to facilitate 
comparison with other exponential approximations 
the exponential approximations were used in the 
integrals 

00 

F(s,r) = / e' 1 rt (l - 5-*- )dt 

s /(1+t ) (10) 
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and 


G(s,r) = / e- irt t(l - 
s 


— i-y ) dt 


( 11 ) 


4 

g(t) - I * k exp (-b k t) (14) 

■> * 1 


!• 


* 


that are related to T and 6 by integration by 
parts (see reference 4 for details of the 
integration by parts). The function to be 
approximated by exponential approximations is 

f(t) = 1 - t//(l + t 2 ) (12) 


The exponential approximations are 


where the a k , b k are tabulated below: 


a! = .101 

a2 = .899 

a 3 = .0047404665 i 

a 4 = -.0047404665 1 

bi = .329 

b 2 = 1.4067 

b 3 = 2.9 + 3.1415926 i 

b 4 = 2.9 - 3.1415926 i (15) 


g(t) = l a. exp (-b k t) (13) 

k=l K K 

This reduces the evaluation of F and G to the 
evaluation of elementary integrals. Equation 
(13) can only be used if t is non-negative. If 
t is negative then g(t) = 2 - g(|t|). 

Replacing f(t) by an exponential approxi- 
mation is not the only way to evaluate F(s,r) 
and G(s,r). Before electing to use this method 
several other techniques were investigated in 
some detail. Direct numerical integration was 
found to be too expensive. A series expansion 
of one or the other integral factor followed by 
a closed form integration is economical; but 
every such expansion has a very limited area of 
application in the s-r plane. To evaluate 
F(s,r) by a combination of an asymptotic series 
for large s (obtained by iterated integration by 
parts) and closed form integration of various 
series expansions of exp (-irt) and t//(l+t 2 ) 
would take at least half a dozen different 
algorithms. This would result in a computer 
subprogram that, while inexpensive to use, would 
be very expensive to develop, certify, and 
maintain. 

For the reasons stated above, exponential 
approximations having the form shown in equation 
(13) are the most commonly used method for 
evaluating the integrals F(s,r) and G(s,r). In 
the next section several currently used approxi- 
mations are discussed. 


This approximation will be referred to as 
W4 (Watkins, 4 term). The error in this 
approximation to f(t), that is e(t) = g(t) - 
f(t), is plotted as short dashes in Figure 1. 

The approximation W4 was derived long before 
the ready availability of high speed computing 
machinery. The coefficients, a(k), and 
exponents, b(k), were generated as follows: 

1. A large scale plot of the function 
f(t) = l-t//(l + t 2 ) was constructed on 
graph paper. 

2. Several plots of a exp(-bt) for different 
values of a and b were drawn on the same 
sheet of paper. Eventually, by trial and 
and error, an a, b were found that appeared 
to give a best fit. 

3. The difference between the function being 
fitted, f(t), and the approximation 

a exp (-bt) was drawn on a large scale plot. 

4. Steps 2 and 3 above were repeated several 
times, using successive difference plots, to 
give an exponential series with several 
terms. 

A more commonly used approximation to f(t) is 
the 11 term exponential polynomial 

11 kbt 

g(t) = l a k e (16) 

k=l K 


A SURVEY OF PREVIOUS EXPONENTIAL APPROXIMATIONS 

There are two reasons for surveying previous 
exponential approximations to f (t ) . One, of 
course, is to point out their deficiencies to 
justify developing a new family of 
approximations. A more important reason is to 
show how they contributed to the development of 
the least squares method described herein and to 
show which features of previous methods have 
been retained. 

The first use of an exponential approximation 
to f(t) to evaluate the Integral F(s,r) was by 
Watkins, Woolston, and Cunningham (reference 1). 
They used the approximation 


due to Laschka (reference 5). The exponent, b, 
and the coefficients, a k , appearing in this 
approximation are tabulated below: 

b = .372 
a x = .24186198 
a 2 = -2.7918027 
a 3 = 24.991079 
a 4 = -111.59196 
a 5 = 271.43549 
a 6 = -305.75288 
a 7 = 41.1836 
a 3 = 545.98537 
ag = -644.78155 
a 10 = 328.72755 

an = -64.279511 (17) 
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This approximation will be referred to as 
Lll. In Lll, instead of II unrelated exponents 
b|<, there is a single exponent, b, and its 
multiples. This means that even though Lll has 
many more terms than W4 the Integral F(s,r) can 
be evaluated almost as fast using Lll as using 
W4 because only a single exponential function 
has to be evaluated. Figure 1 (long dashes) 
also shows that Lll is more accurate than W4. 

The computer code to approximate the integral 
F(s,r) is somewhat easier to write using the 
approximation Lll than it is using the 
approximation W4. This is probably why the 
approximation Lll is so widely used. 

Reference 5 does not state how the exponent, 
b, qr the coefficients, a^, in Lll were 
computed. Inspection of Figure 1 leads one to 
suspect that they were computed in essentially 
the same manner as the coefficients and 
exponents of W4. There are 12 free parameters 
in this approximation, 1 exponent and 11 
coefficients. If they had been computed by some 
sort of systematic optimization procedures such 
as least squares or minimax then the number of 
nodes (points where e(t) = 0) on the error plot, 
Figure 1, would approximately equal the number 
of free parameters. 

The approximation Dll has been In use for 
several years in a computer program based on 
reference 1 but the coefficients have not been 
published. It has the same form as Lll, namely 
equation (16). However the coefficients and 
exponent were computed to minimize the 
integrated square error 

E = f w(t) (g(t) - f(t)) 2 dt (18) 

o 

using the weight function w(t) = t-1/2. The 
reason for this choice will be discussed later. 
The exponent and coefficients for Dll are 
tabulated below: 


b = .1539 
aj = .0545 
a 2 = -.8049 
a 3 = 6.3712 
a 4 = -10.8582 
a 5 = -130.6866 
a 6 = 915.2209 
■a 7 = -2793.8215 
a 8 = 4800.8897 
ag = -4792.6556 
a 10 = 2596.7163 
an = -590.2260 

The least squares procedure used to generate 
Dll was very easy to implement because there is 
only one non-linear parameter, b. The form of 
the approximation Lll was chosen to minimize the 
execution time of the computed approximation but 
it also makes the coefficients in the 
approximation easier to compute. 

Figure 1 (solid curve) indicates that Dll is 
much more accurate than Lll. Since it has the 
same form as Lll a program that uses Lll can be 
converted to one using Dll by only changing one 
or two data statements. This similiarity in 


form is the only reason for tabulating the Dll 
coefficients because It will be shown later that 
much more accurate approximations can be ob- 
tained by changing the expression for the 
exponents b^ to something other than k*b. 

Before considering this, however, it's instruc- 
tive to see what happens to an approximation 
such as Dll when the number of terms, n is 
increased. 

Consider the approximation 


g(t) = l a. e‘ kbt (20) 

k=l 

where the exponent b and the coefficient a^ 
are chosen to minimize the integral E of 
equation (18). Figure 2 shows the e(t) vs. t 
plots for three of these approximations. For 
the three approximations shown the maximum error 
decreases with n as expected. The ragged 
behavior of the 24 term approximation for t 
near zero is due to cancellation errors. The 
coefficients a^ are very large while both the 
function being approximated f(t) and the approx- 
imation functions exp(-kbt) are near unity. For 
n slightly greater than 24 this cancellation 
due to large coefficients becomes the main 
source of error in these approximations and no 
further gain In accuracy is attainable using 
b|< = k*b unless the computer word length is 
increased. 

The solution to this cancellation error 
problem is suggested by an exponential 
approxiation due to Jordan (ref. 6). The e(t) 
vs. t plot for this approxiation, called J10, is 
shown as the solid curve in Figure 3. 

This approximation is of the form 

g(t) = a exp(-b t) + l a k (21) 

o o k=1 x 

with exponents and coefficients tabulated below: 
bo = 3 

a 0 = .7048426 
b = .0625 
a! = .002907843 
a 2 = .002591528 
a 3 = .02667074 
a 4 = .070971 
a 5 = .347837 
ag = .5556069 
a 7 = -.776979 
a 8 = .07004561 

. a 9 = -.004557519 (22) 

Although the exponents bk in this 
approximation are somewhat more complicated than 
in Lll or Dll, the sum in equation (21), when 
substituted into the integral F(s,r), can still 
be evaluated without repeated calls to the 
exponentiation subprogram. This simplification 
occurs because exp(-2*bt) is equal to 
(exp (-2(k-l).b*t)2 so that the exponential 
functions can be computed recursively as in Lll. 

The approximation J10 is much more accurate 
than Lll and is somewhat more accurate than 


* 


# 

% 
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Dll. It would take a 15 term approximation 
similar to Dll to attain the accuracy of J10. 
However it isn't the accuracy of J10 that is 
important. It is the fact that the 
coefficients, a^, in J10 are all bounded by 
unity so cancellation errors are very small. 

Reference 6, like reference 1 and 5, does not 
state how the exponents and coefficients were 
computed. The coefficients in J10 were probably 
computed in a manner similar to that used for 
the W4 coefficients. This is indicated by the 
fact that the number of nodes shown in figure 3 
(eight) is much less than the number of free 
parameters (twelve). By contrast, the 
approximations D8.1 and D12.1, developed using 
the methods of the next section, and also 
plotted in figure 3, have the same number of 
nodes as free parameters. 


EXPONENTIAL APPROXIMATIONS OF ARBITRARY ACCURACY 

In this section a family of new exponential 
approximations of the form 

g(t) = l a k exp (2 k /">bt) (23) 

k=l K 

will be investigated. For each of these 
approximation the exponent b and the 
coefficients ak were computed to minimize the 
integral E of equation (18). The coefficients 
ak enter into the error minimization process 
linearly and the exponent b enters 
nonlinearly. As a consequence, for each value 
of n and m there are several values of b for 
which E is a relative minimum. Of these 
relative minima the value of b chosen was not 
the value which minimized E in the absolute 
sense but rather, the values that minimized the 
maximum values of |e(t ) j = |g(t) - f(t)|. 

Exponents and coefficients for four of these 
approximations, designated D8.1, D12.1 , D24.2, 
and D72.3 (Dn.m) are tabulated in TABLE I. The 
error e(t) in D8.1 and D12.1 is also plotted in 
Figure 3. A summary of some of the properties 
of these approximations is listed in TABLE II. 
Properties of the approximations discussed 
previously are also included in TABLE II. 

Approximations D8.1 and D12.1 can be used to 
rapidly evaluate the integral F(s,r). The 
execution times (TABLE II) are comparable to 
those of approximations, W4, LI 1/ Dll , and J10. 
However D8.1 and D12.1 are much more accurate. 
The author knows of no study to determine how 
accurately the integral F(s,r) has to be 
computed for various aeroelastic applications. 
Approximation D24.2 permits computing the 
integral F(s,r) to about 3.5 more significant 
figures than D12.1. This is about full register 
accuracy on a short word length computer. 
Approximation D72.2 is probably more accurate 
than is needed for any engineering purpose. 

These coefficients are included because they 
were used to assess the error of the other 
approximations. 


The form of the expression chosen for the 
exponents, namely 

bk = 2 k/m b . (24) 

was dictated to a large extent by the approx- 
imations discussed in the previous section. In 
Lll the form k-b was introduced to reduce the 
number of exponential function evaluations that 
occured when evaluating F(s,r). However this 
choice also reduced the number of nonlinear 
parameters in the least squares error mini- 
mization process used to generate Dll. This 
made it possible to automate the process. The 
choice, bk = 2 k *b, suggested by J10, is 
almost as good as that of Dll for convenience in 
computing the integral F(s,r) and is better for 
the least squares process. This is because the 
matrix of the normal equations has a much lower 
condition number when 2 k *b is used than it has 
when k*b is Used. 

The choice, bk = 2 k *b, corresponding to 
m = 1 in equation (23), was used to generate the 
approximation D8.1 and D12.1. However when an 
attempt was made to generate D24.1 this way the 
result was only slightly more accurate than 
D12.1. The exponent, b enters into the least 
squares process in a nonlinear manner so there 
is more than one value of b for which E = E(b) 
is a relative minimum. When n is 12 there are 
eight of these relative minima and the lowest 
value of E is 1.56E-9. Increasing n to 24 
increases the number of relative minima to 
twenty but only reduces the lowest value of E to 
9.07E-1Q. In fact, when n is. 24 there are ten 
relative minima with integrated square errors, 

E, between 9.07E-10 and 9.9E-10. As n is 
increases a value of n is reached for which 
increasing n increases the number of relative 
minima without decreasing the value of E at the 
lowest relative minimum. This happens because 
when n is increased from 12 to 24 the value of b 
was rescaled so that effectively the additional 
bk were all inserted ahead of bi, where 
their effect was to reduce the te (t ) j for t near 
zero, and they were also inserted after b^g, 
where their effect was to reduce |e(t ) j for t 
very large. None of them were interpolated 
between the bk to reduce le ( t ) | in the middle 
of its range. Changing b|| from 2 k *b to 
2 k /2.b has the desired effect of reducing 
J e ( t ) ( in the middle of its range. When n is 
set to 24 and m is set to 2 only seven relative 
minima are obtained and the value of E at the 
lowest is 1.78E-12, a significant improvement 
over m = 1. One might think that since m = 2 is 
an improvement when n = 24, that m = 3 would be 
even better. When m was set to 3 the number of 
relative minima was reduced to four but the 
lowest value of E was increased to 1.9E-10. 

Also for n = 24 and m = 3 the coefficients ak 
were very large so some cancelation errors can 
occur. For most values of n there is a single 
optimum value of m. If there are two optimum 
values of m for a particular value of n then the 
lower should be chosen because evaluating 
F(s,r) requires m exponential function 
evaluations. 

Exponential approximations of the form used 
in equation (23) can be used to generate 
approximations to f(t) that have arbitrary 
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accuracy. All that has to be done is to keep 
increasing the number of terms n and whenever 
this fails to increase the accuracy, then 
increase the separation parameter m. There are 
no cancelation problems due to the a^ becoming 
extremely large as happens when b|< = k*b. The 
normal equations of the least squares process 
have matrices whose condition number increases 
with n but the rate of increase is very slow. 

DISCUSSION 

The following recommendations regarding the 
use of exponential approxiamtions are based on 
the information in table II: 

1. If execution time is very important use 
approximation D8.1. 

2. Any programs for which execution time is not 
critically important that currently use 
approximations W4 or Lll should be changed to 
use approximation D12.1. This will reduce the 
error by two order of magnitude at the expense 
of a negligible increase in execution time. In 
fact, it was found when converting a particular 
program from Lll to D12.1 that the execution 
time was decreased. This was because during the 
recoding some compiler generated complex 
arithmetic was replaced by separate evaluation 
of real and imaginary parts. 

3. New programs should use approximation D24.2. 
It gives full word length accuracy on short word 
length computers and takes less than twice as 
long to execute as 012.1. 

4. Programs that currently use approximation 
J10 should only be converted to D12.1 if a 
considerable amount of future use of the program 
is anticipated. Approximation D12.1 is no 
faster than J10 and Is only an order of 
magnitude more accurate. 

In non-planar flow two incomplete modified 
Struve functions are required. They both use 
the same exponential approximation but the 
non-planar term has an extra variable of 
integration factor. This acts as a weight 
function on the error. As a consequence the 
exponential approximations described in this 
report introduce more error into the non-planar 
part of the kernel than into the planar part. 
This was not investigated in any detail. It is 
felt that it is not very important because all 
successful non-planar programs, such as the 
various doublet lattice programs, use very large 
panels so the kernel is either planar or is 
attenuated by distance. 

The integrals F(s,r) and G(s,r) are integrated 
by merely substituting g(t), defined by equation 
(23), for f(t), in equations (10) or (11) and 
then integrating termwise. Thus, if s 0 

n 

F(s.r) - e' 1rS l a k e k /(2 k/ Vir) (25) 

k=l K K 

G(s,r) « sF(s,r) - e' irs l a,e. /(2 k/,m b+ir) 2 (26) 

k=l K K 


where 

e k = exp(-2 k/m bs) (27) 

The exponential factors In equations (25) and 
(26) should be computed from equation (27) only 
if k < m. If k is greater than m then e^ 
shoulcf be computed from the recursion 


Occasionally it is necessary to compute F(s,r) 
and G(s,r) for complex r. This occurs when the 
frequency u> in equation (1) has an imaginary 
component. If Imfr} < 0 then the integrals 
(10) and (11) exist and the sums (25) and (26) 
are valid approximations. A minor problem is 
that the error estimates given by the third 
column of table II may be too low because they 
are based on a search of the s-r plane for real 
s and r. If Im { r } > 0 then the integrals (10) 
and (11) do not exist. This is because the 
oscillations are damped and the integrals (10) 
and (11) usually represent integration over the 
subsonic wake. Even if Im(r) > 0 the sums (25) and 
(26) exist. They are approximation to Abel sums 
of the divergent integrals (10) and (11). The 
Abel sum of an exponentially diverging Integral 
is obtained by inserting a factor exp(-et2) 
into the integrand and taking the limit as e 
approaches zero. If it can be shown on physical 
grounds that the integrals should exist then 
they are defined by their Abel sums. 

CONCLUDING REMARKS 

Several different methods for evaluating the 
incomplete modified Struve functions that occur 
in unsteady aerodynamics were investigated. Of 
the methods investigated it was found that an 
exponential approximation to the algebraic part 
of the defining integral followed by a closed 
form integration furnished the best compromise 
among accuracy, execution speed, and code 
maintainability. If the exponent spacing of the 
exponential approximation is either an 
arithmetic sequence, as suggested by Laschka, or 
a geometric sequence, as suggested by Jordan, 
then the number of exponential function 
evaluations required to evaluate the incomplete 
Struve functions is reduced. This same sort of 
spacing also simplifies the evaluation of the 
coefficients of the approximation if least 
squares are used. Of the two spacing 
sequences the geometric is much better. For 
either spacing sequence least squares provides 
an easily automated way of computing the 
coefficients, a^ and the sequence multiplier 
b. Approximations with coefficients computed by 
least squares are more accurate than those of 
the same order with coefficients computed by 
manual curve fitting. 

In the paper a very limited class of 
exponential approximations was Investigated, 
namely those with arithmetic and geometric 
spacing exponent sequences and with coefficients 
computed by least squares. It is possible that 
other sequences and schemes for computing the 
coefficients would give more accuracy. In 
particular using a minimax algorithm to generate 
the coefficients and exponent multiplier would 


6 



certainly give a small Increase In accuracy. A 
program to Implement a minimax algorithm is much 
harder to write and is much more expensive to 
execute than a least squares program. 

Furthermore the Increase In accuracy would only 
be very slight. However if a minimax program 
were available it could be used to generate 
exponential approximations over subintervals of 
the positive real axis with no greater effort 
than is required to fit the whole axis. This is 
not true for least squares because of the 
integrals that have to be computed (the ones 
that contained quarter order Bessel functions). 
Splitting the interval into subintervals can 
lead to a large increase in accuracy with no 
increase in computing time. The only penalty is 
in program maintenance and development cost. 
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TABLE I 


Exponent multiplier, b and coefficients, 
ak, for several approximations to 
l-t//(l+t 2 ) of the form indicated by equation 
(23). 



Approximation D8.1 


n = 8 

m = 1 

b = 

.035003907466 

a l = 
a 3 = 
a 5 = 
a 7 = 

.004329519485 

.033195062769 

.376739860841 

-.380262739620 

a 2 = 
a 4 = 
a 6 = 
a 8 = 

.001601370746 

.098682301170 

.822464185014 

.043400039240 


Approximation D12.1 


n = 

12 

m = 1 

b 

= 

.009054814793 

a l 


.000319759140 

a 2 

— 

-.000055461471 

a 3 


.002726074362 

a 4 

= 

.005749551566 

a 5 

= 

.031455895072 

a 6 

= 

.106031126212 

a 7 


.406838011567 

a 8 

ss 

.798112357155 

ag 

= 

-.417749229098 

a 10 

s 

.077480713894 

a ll 


-.012677284771 

a 12 

=t 

.001787032960 


Approximation D24.2 


n = 

24 

m = 2 

b 

s 

.005209230865 

a l 

— 

.000305311497 

a 2 

— 

-.001412280807 

a 3 

= 

.003845227615 

a 4 


-.007196572664 

a 5 

= 

.011385147609 

a 6 

= 

-.014763498650 

a 7 

= 

.018969114027 

a 8 

= 

-.019842326360 

a 9 


.025618710871 

a 10 

S 

-.020313397232 

a ll 

= 

.036575115249 

a 12 

= 

-.010202806435 

a 13 

= 

.069407344423 

a 14 


.037308217964 

a 15 

s 

.177803740980 

a 16 

= 

.198282197469 

a 17 

= 

.433959048197 

a 18 

= 

.354218469431 

a 19 

= 

.104676453558 

a 20 

3 

-.715978991168 

a 21 

= 

.407542943867 

a 22 

3 

-.104393578248 

a 23 

55 

.015398943987 

a 24 

3 

-.001192670868 


Approximation D72.3 


n = 

72 

m = 3 

b 

= 

.000065986269 

a l 

— 

.000000487572 

a 2 

= 

-.000003844799 

a 3 

3 

.000015710073 

a 4 

a 

-.000044143564 

a 5 

= 

.000096360019 

a 6 

=s 

-.000174937155 

a 7 

3 

.000276395746 

a 8 

= 

-.000392188471 

a 9 

= 

.000511902333 

a 10 

a 

-.000625480027 

a ll 

= 

.000725938341 

a 12 

= 

-.000808476891 

a 13 

3 

.000872553959 

a 14 

= 

-.000917456689 

a 15 

= 

.000947455433 

a 16 

= 

-.000961741082 

a 17 

= 

.000968885391 

a 18 

= 

-.000962841735 

a 19 

= 

.000960418999 

a 20 

= 

-.000941494817 

a 21 

= 

.000943838490 

a 22 

* 

-.000912640747 

a 23 

= 

.000939494732 

a 24 

=s 

-.000883200418 

a 25 

= 

.000971868480 

a 26 

= 

-.000847331705 

a 27 

3 

.001082186979 

a 28 

= 

-.000771912330 

a 29 

= 

.001357091273 

a 30 

= 

-.000556797879 

a 31 

3 

.001997471929 

a 32 

= 

.000067960260 

a 33 

= 

.003487618691 

a 34 

= 

.001801105035 

a 35 


.007010443761 

a 36 


.006406650025 

a 37 

= 

.015440810290 

a 38 

= 

.018199043007 

a 39 

= 

.035536908427 

a 40 

= 

.047032188464 

a 41 

= 

.081594062558 

a 42 

= 

.111356164158 

a 43 

3 

.174138343702 

a 44 

= 

.222492227365 

a 45 

= 

.288676153073 

a 46 

= 

.263088797407 

a 47 

= 

.143795607125 

a 48 

= 

-.194655408459 

a 49 

3 - 

-.414538285804 

a 50 

a 

-.093514761246 

a 51 

= 

.562053915950 

a 52 

= 

-.395454683729 
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a 53 = .132558644137 
a 55 = -.022964836837 
a 57 = -.024548039318 
a 59 = -.017250853224 
ag} = -.011265419745 
a 6 3 = -.006583208288 
a 65 = -.003149336650 
a 67 = -.001071671331 
a 69 = -.000203826307 
a 71 » -.000012950446 


a 54 = -.010481139795 
a 56 = .027218864137 
a 58 = .020804555835 
a 60 = .014078375164 
a 6 2 = .008772771052 
a64 = .004702560548 
a 6 6 = .001939706215 
368 = .000513691017 
a 7 o = .000062322287 
a 7 2 = .000001360075 


TA8LE II 

Properties of various approximations 
to 1 - t//(l + t2) . 


In the least squares procedure the coef- 
ficients ak and the exponent b in equation 
(A2) are chosen to minimize the Integral 

E = f w(t) [g(t ) - f(t)] 2 dt (A3) . 

o 

There is no reason to suppose that this is an 
optimum choice for the a^ and b. In fact, 
there is some evidence the minimizing the 
maximum error j g ( t ) - f(t)| would be a better 
choice. However, minimizing E certainly gives a 
good set of coefficients and exponent and the 
process is very easy to automate. If the weight 
function w(t) in equation (A3) is set to 1 then 
the absolute value of the error function 



max )e(t)l 

max | E ( s , r ) J 

time 

W4 

1.6E-3 

2.3E-2 

355 

Lll 

1.3E-3 

1.8E-2 

382 

Dll 

3.5E-4 

2.6E-3 

382 

J10 

1.3E-4 

1.7E-3 

407 

D8.1 

1.6E-4 

1.1E-3 

293 

D 12.1 

2.5E-5 

1.9E-4 

411 

D24.2 

3.5E-7 

2 . IE -6 

768 

D72.3 

3.0E-10 


2250 


The first column of the table is the 
approximation designator. 


|e(t)| «= | g ( t ) - f(t ) | (A4) 

takes on its maximum values near the ends of the 
approximation interval. The weight function 

w(t) = (t-a)' 1/2 (b-t)' 1/2 (A5) 

is frequently used to make a least squares 
approximation over a finite interval (a,b) 
closely resemble a minimax approximation. In 
equation (A3) b = « and a = 0 so the (b-t ) -1/2 
factor is undefined and the (t -a ) - 1/2 factor 
i s 


The second column is the maximum error in the 
integrand approximation. It is very easy to 
compute but not too closely related to the error 
in the integral. 

The third column is the maximim error in the 
integral for the planar component only. It is 
extremely expensive to compute. 

The fourth column is the computer execution 
time for a single evaluation of the planar 
integral with positive lower limit and infinite 
upper limit. It is expressed as a multiple of 
the computer's add time. 


APPENDIX 

THE LEAST SQUARES PROCEDURE FOR COMPUTING 
COEFFICIENTS AND EXPONENTS 


w(t) = t" 1/2 ( A6) 

Use of this weight function also makes some 
required integrals easier to evaluate. 

The coefficients a^ and the exponent 
multiplier b are computed by setting aE/aa^ to 
zero for any value of b and by setting dE/db to 
zero where the ak are functions of b defined 
by the condition 3E/ aa^ = 0. 

l.t 7 H t ‘ » <«> 

Then 

/ t_1/ * [9(t) - f (t) ] exp(-bp A t) dt = 0 (A8) 

That is 


The function 

f(t) = 1 


t 

/(I + t 2 ) 


(Al) 


is to be approximated by an exponential sum 


n 

g(t) = l a k exp(-bp.t) (A2) 

k=l K K 


The number of terms, n, is preassigned as are 
the exponent coefficients Pk. Although the 
only pk arrays for which calculations have 
been performed are Pk = k and pk = 2 k / m 
the least squares procedure is valid for any 
preassigned pk array. The only restriction is 
that all the pk be distinct. 


I a / t" 1/2 exp[-b(p + p.)t] dt 
k=l K o 1 K 

0 o 

= / r 1/Z exp(-bpt) f(t) dt - (A9) 
o * 

The Laplace transform of fl/2 i s 

/ e' pt f 1/2 dt = /„ p-1/2 (A10) 

o 

so equation (A9) can be written 
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2 =«k*k ’ " t « A,1 » 

k=o 


where 

c .k • o, * >v >‘ 1/2 < #12 ) 


where a'k = da^/db is obtained by 
differentiating equation (All) and solving 


n da k . bp 
Skiff) = “27(1^) H 2* 

and where 



. bp 

/(b/»)H ( — |) 
(A21) 


d £ = /(b/n) H(bp £ /2) (A13) 


H* (y) = - jy H(y) + /(Try/2) 




and 

H(y) = /(^y) - f e' 2yt t 1/Z (l + t 2 )- 1/Z dt (A14) 

A closed form expression for H(y) is obtained by 
differentiating the Laplace transform pair 
4.3(32) of ref. 7 and then exploiting some 
relationships between parabolic cylinder 
functions and Bessel functions, equations 19.3.7 
and 19.15.9 of ref. 8. The expression is 

H(y) = *(fy) 1/2 [J 1/4 (y) J. 3/4 (y) 

■ J _i/4^ y ) J 3/4^ y ^ + /z J l/4^) J 3/4^ y ^ 

(A15) 

After the ak are compute then E can be 
evaluated as a function of b. 

E - / f 1/z (g-f) g dt - / f 1/2 (g-f) f dt 
o o (A16) 

The first integral in (A16) is zero because of 
equation (A8) so 


+ J 3/4^ ‘ J l/4^ + 1/2 f J -3/4^ y ^ J 3/4^ y ^ 


♦ J 1/4 (y) J. 1/4 (y)]} («2) 

A detailed description of the algorithms used to 
implement the above equations follows: 

I. Since E 0 is independent of n and b it 
only has to be computed once. The quantity 

r ( 1/ 4) appearing in equation ( A19) was computed 
by first using equations 6.1.48 of ref. 8 to 
compute r( 145+1/4) and then using the functional 
equation for r(z), equation 6.1.15 of ref. 8, to 
compute r( 1/4) . 

II. The number of terms, n is read into the 
program. The elements of the matrix C =[ c zk]» 
which are Independent of b, are computed from 
equation (A12). Since C is symmetric only the 
diagonal and upper triangle of C are stored. 

The C matrix, which is positive definite, is 
factored without pivoting 


C = IDU ( A23) 


where 


E = E o - I a k H(bp k /2) 

k=l 


E 0 = / t‘ 1/2 [f(t)] 2 dt 
0 


(A17) 


(A18) 


This can also be integrated in closed form 


71 r 8/ ( 2lT ) 

* 71 (r(l/4)) 2 


-1] 


1.16741087 

(A19) 


The exponent b is obtained by solving dE/db = 0 


§ = - I [a' H(bp k /2) + a k (p k /2) H*(bp k /2) 

(A20) 


where L is a unit («kk = 1) lower triangular 
matrix, D is a diagonal matrix, and U is a unit 
upper triangular matrix. Since L is the 
transpose of U only D and U are stored. They 
are overwritten onto C in memory. The algorithm 
dependent spectral norm condition number of C, 
max j d^k | /mi n |dkk|» is also computed. 

III. For each value of b both E(b) and 
E '(b) = dE/db are computed from equations (A17) 
and (A20) . The a|< and a'k are obtained by 
solving equations (All) and (A20) using the LDU 
factorization described in step II above. The 
quantities H(y) and H'(y) for y = b*pk, k=l to 
n, are computed from the four Bessel functions 
J-3/4, J-l/4> J l/4> and J 3/4- These 
Bessel functions are computed using the 
Miller-Abramowitz algorithm described in example 
1 of section 9.12 of ref. 7. The example is for 
integer order and uses equation 9.1.46 for 
normalization. Since Ji/ 4 , etc. are of 
non-integer order equation 9.1.87 of ref. 7 has 
to be used for normalization. The highest 
argument gamma function in eq. 9.1.87 is 
computed as in Step I and the rest of the gamma 
functions are computed recursively from eq. 
6.1.15. 


IV. Finding a value of b for which E(b) is a 
minimum is performed in two steps. 
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A. A plot of logigE vs. logjob is 
constructed in order to determine the number of 
minima and their approximate locations. 

B. Each approximate minimum is identified by 
three consequtive b values, bj, bp, and 

b3, for which the central value of E(b), Ep 
is lower than either of the two adjacent values, 
E i , and E3. These three values of b and the 
associated values of E'(b), E'j, E ' 2 » and 
E ' 3 are used as starters for Muller Iteration 
which is used to solve the equation E'(b) = 0 . 


V. Step IV above determines many values of b 
for which E(b) is a relative minimum. For each 
of these values of b the maximum value of |e(t)l 
was obtained from a plot of e(t) vs t. The plot 
for which maxle(t)| was a minimum gave the 
optimum value 'of b! 

The program to implement the least squares 
process was originally written for p|< = k. 

For this choice of Pk array the C-matrix 
(C£k = (4 + k)-l/ 2 ) is extremely ill 
conditioned so all calculations were performed 
in double precison. Double precision is not 
needed for pk = 2^/ (n but was used because 
the program had already been written. 


10 



Error e(t) = g(t) - f(t) 



0 2 5 1 0 20 50 o 

t 

Figure 1. - Comparison of error in the exponential approximations W4, Lll, and 
Dll to the function l-t//(l+t 2 ) . The nonlinear t-scale on the plot results 
from mapping the (0. -) t-interval into a (-1,1) u-lnterval by the bilinear 
transformation u = (t+10)/(t-10). 
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Figure 2. - Errors In least squares approximations for three values of n 
using Laschka's form for the exponents. Notice that the vertical scale Is 
different from that of figure 1. 
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